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Abstract 

Using molecular dynamics simulations we investigate the structure of a system of particles 
interacting through a continuous core-softened interparticle potential. We found for the trans- 
lational order parameter, t, a local maximum at a density pt-max and a local minimum at 
Pt-min > Pt-max- Between pt-max and pt-min, the t parameter anomalously decreases upon in- 
creasing pressure. For the orientational order parameter, Qq, was observed a maximum at a 
density pt-max < PQmax < Pt-min- For densities between pQ max and p t -min, both the translational 
(t) and orientational (Qe) order parameters have anomalous behavior. We know that this system 
also exhibits density and diffusion anomaly. We found that the region in the pressure-temperature 
phase-diagram of the structural anomaly englobes the region of the diffusion anomaly that is larger 
than the region limited by the temperature of maximum density. This cascade of anomalies (struc- 
tural, dynamic and thermodynamic) for our model has the same hierarchy of that one observed 
for the SPC/E water. 



PACS numbers: 64.70.Pf, 82.70.Dd, 83.10.Rs, 61.20.Ja 



I. INTRODUCTION 



Water is the most important substance for life: It cools, carries, stabilizes, reacts, lu- 
bricates, dilutes, and much more. Despite of this, many of its characteristics are not well 
understood. While most liquids contract upon cooling, water expands below T = A°C at 
ambient pressure^. This is known as the density anomaly of water. Heating the water from 
T = 0°C up to T = 4°C a competition between open low density and a closed high density 
structure takes place. The gain of thermal energy breaks a considerable number of hydrogen 
bonds what leads the open low density structure to become unstable in relation to the closed 
high density structure. So, the system contracts. 

Density anomaly is not the only one, far from it, the literature reports forty-one anomalies 
for wateA Not only the thermodynamics of water is anomalous, but also its dynamics. 
Commonly the materials diffusivity decreases with increasing pressure. Liquid water has an 
opposite behavior in a large region of the phase diagra m 3 i 4 i 5 i 6 i 7 i 8 i 9 i 10 i 1:L . Increase in pressure 
disturbs the structure by inclusion of interstitial molecules that share an hydrogen bond 
with another one. As a result, the bond is weakened and the molecule is free to move. The 
shared bond breaks and the molecule by means of a small rotation, connects to another 
molecule enabling the translational diffusion 3 -. 

Water is not an isolated case. There are other examples of tetrahedrally bonded molecular 
liquids such as silica and silico n 12 ! 13 that exhibit thermodynamic and dynamic anomalies. 
Thermodynamic anomalies were also found in liquid metals^ and graphite^. Unfortunately, 
a closed theory giving the relation between the form of the interaction potential and the 
presence of the anomalies is still missing. 

It is reasonable to think that the structure and anomalies are deeply related. Establishing 
the connection between structure and the thermodynamic and dynamic behavior of water is 
a fundamental step towards understanding the source of the anomalies. At this point a ques- 
tion emerges: how can we define (measure) structure in liquids? Errington and Debenedetti^ 
proposed two simple metrics: a translational order parameter, t, that measures the ten- 
dency of pairs of molecules to adopt preferential separations, and the orientational order 
parameter^£, q, quantifying the extend to which a molecule and its four nearest neighbors 
assume a tetrahedral arrangement. For other crystal configurations one may use the orienta- 
tional order parameter introduced by Steinhardt et a/.—, Qq, which depends on the number 
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of nearest neighbors taken into account for each molecule. For a completely uncorrelated 
system (ideal gas) both t and q must to be zero and Qq, is equal to one over the square of 
the number of neighbors. For a crystal, t, q and Q 6 are large. Torquato el alM introduced a 
systematic way to study the structural order in liquids mapping state points into the t — q 
plane. They refer to it as an order map. Errington and Debenedetti used the order map to 
investigate structural order in simple point charge/extended (SPC/E) water 8 . 

For normal liquids, t and q increase upon compression, because the system tends to be 
more structured. It was found that in SPC/E water both t and q decrease upon compression 
in a certain region of the pressure-temperature (P-T) phase diagram 8 . This region is referred 
as the region of structural anomalies. Errington and Debenedetti showed also that, inside 
the structurally anomalous region, all the paths formed by the (t, q) points collapse into a 
single line. This means that the translational order parameter, t, and the orientational order 
parameter, q, are coupled. Outside the structurally anomalous region the states points in 
the order map define a two-dimensional region, meaning that the parameters t and q are 
independent. 

Performing molecular dynamics simulations, Errington and Debenedetti and Netz et a/.— 
showed that in SPC/E water the thermodynamic and dynamic anomalies form nested domes 
in the P-T phase- diagram, where the diffusion anomaly lies outside the density anomaly. 
Additionally, Errington and Debenedetti showed that the structurally anomalous region 
englobes the diffusion and density anomalies regions. 

Several models of water for computer simulations have been proposed^, with three, four 
or five localized partial charges, some of them having Lennard- Jones interaction centers in 
the oxygens and hydrogens, others only in the oxygens. A considerable number of these ap- 
proaches reproduce many anomalies present in liquid water. However, these models are 
complicated, what makes difficult to understand the physics behind the anomalies. In 
this sense, isotropic models are the simplest framework to understand the physics of liq- 
uid state anomalies. Moreover, the use of an effective potential is particularly suitable 
for extending our conclusions for more complex fluids. From the desire of constructing 
a simple two-body isotropic potential, capable of describing water-like anomalies, a num- 
ber of models in which single component systems of particles interact via core-softened 
(CS) potentials^ have been proposed. They possess a repulsive core that exhibits a re- 
gion of softening where the slope changes dramatically. This region can be a shoulder or a 
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ramp 



22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38.39,40,41,42,43 



In the shoulder case, the potential consists of a hard core, a repulsive shoulder and, 
in some cases, an attractive square we n 22 i 2 3i 2 VV 6 i 27 i 28 i 2 V°i 31 i 32 i 33 i 34 i 37 . rp^g potential has a 
change in the slope at short-ranged distances. In two dimensions, such potentials have ther- 
modynamic and diffusion anomalies. In three dimensions, no dynamic and thermodynamic 
anomalies were reporte d 24 > 25 i 26 i 27 i 31 i 32 . 

In the ramp case, the interaction potential has two competing equilibrium distances, de- 
fined by a repulsive ramp 36 ' 37 ' 38 ' 40 ' 41 ' 42 ' 43 . In some cases an attractive part is include d 36 ' 3 ^ 42 . 
In two dimensions, there are thermodynamic anomalies in such potentials. In three dimen- 
sions, these potentials exhibit not only thermodynamic anomalies, but also dynamic and 
structural anomalie o 38 i 40 i 41 i 43 . 

Notwithstanding the progresses described above, a model in which both the potential and 
the force are continuous functions and that exhibits all the thermodynamic and dynamic 
anomalies like the ones present in water is still missing. In this paper, we check if a ramp-like 
potential previously studied by us^ 3 - has not only density and diffusion anomalies, but also 
structural anomalies. We will verify if the region in the pressure-temperature phase-diagram 
of thermodynamic and dynamic anomalies are inside the region of structural anomalies as in 
SPC/E water—. The hierarchy between the anomalies in such simple model is an important 
step in order to understand the mechanism of the anomalies. 

The reminder of this paper goes as follows. In sec. II the model is introduced; in sec. 
Ill the methods for calculating structural order in liquids are presented. Results for the 
structural anomalies and the order map obtained from molecular dynamics simulations are 
shown in sec. IV. Conclusions about the relation between the locus of the thermodynamic, 
dynamic and structural anomalies and about the order map are presented in sec. V. 

II. THE MODEL 

The model we study consists of a system of N particles of diameter a interacting through 
an isotropic effective potential given by 
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where U* (r) = U(r)/ e. The first term of Eq. ((J) is a Lennard- Jones potential of well depth 
e and the second term is a Gaussian centered on radius r = with height a and width c. 
Depending on the choice of the values of a, r and c, this potential assumes several shapes 
ranging from a deep double well potentia l 44 ' 45 ' 46 to a repulsive shoulder—. 

Recently, using molecular dynamics simulations and integral equations theory, we have 
studied the potential Eq. (JTJ) setting a = 5, r /a = 0.7 and c = 1 (see Fig. [TJ) 43 <. Here, we 
use the same parameters (as in the previous 4 ^). It is interesting to note that close to the core 
(r/a ~ 1) this potential experiences an unusual change of slope, weakening the repulsive 
force between the particles. 
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FIG. 1: Interaction potential eq. (JJ) with parameters a = 5, r^ja = 0.7 and c = 1, in reduced 
units. The inset shows a zoom in the very small attractive part of the potential. 



III. THE METHODS 



A. Translational order parameter 



The translational order parameter of a system of particles of density p = N/V, where N 
is the number of particles and V is the volume of the system, is defined a o^ 1 ? 116 , 

t= [ C \g(t)-l\dt, (2) 



where £ = rp 1 / 3 is the interparticle distance r divided by the mean separation between pairs 
of particles p" 1 / 3 . g(£) is the radial distribution function, where g is proportional to the 
probability to find a particle at a distance £ to another particle placed at the origin. £ c is 
a cut-off distance. In this work, we use 4 * 7 £ c = p 1 ^ 3 L/2, where L = V 1 ^ 3 . For a completely 
uncorrelated system (ideal gas) g = 1 and t vanishes. In a crystal, a translational long-order 
(<7 7^ 1) persists over long distances making t large. 

B. Orientational order parameter 

For the orientational order parameter introduced by Steinhardt el. alM, we follow the 
strategy introduced by Yan el. alM. We define k vectors, r^, connecting the particle % 
with its k nearest neighbors j. Each vector is a "bond". A polar (4>ij) and azimuthal 
(Oij) angles with reference to an arbitrary axis may be associated to each bond and 
the spherical harmonics Yi m (9ij,<pij) may be calculated. After computing the average of 
Yim(8ij j <f>ij) over the k bonds, namely, 

k 

(t) = lEU%,^), (3) 

J'=l 

one can evaluate the orientational order paramete r 8 i 13 i 16 i 19 i 48 i 49 associated to each particle i, 

1/2 

(4) 

For characterize the local ordep^ of the system was use d 1 

1 - 

i=l 

that is the mean value of Q l 6 over all particles of the system. The Q§ parameter assumes its 
maximum value for a perfect crystal and decreases as the system becomes less structured. For 
a completely uncorrelated system (ideal gas) Q l § =1/ Vk. For a crystal, the Q e value depends 
on the specific crystalline arrangement and the number of neighbors taken into account. For 
example, for the face centered cubic (fee) with its twelve first neighbors (k — 12), we have 
Q{ CC = 0.574. For a body centered cubic (bec), which have only eight nearest neighbors, 
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Qq = 0.628. Note that if we include not eight, but fourteen neighbors for calculating 
Q b 6 cc ~ k , we have Q b 6 cc ~ u = 0.510. 

For the potential given by the Eq. (JTj), the expected crystalline configuration at the 
ground state for low densities is the hexagonal close packing (hep), which have twelve first 
neighbors (see sec. II VI for more details). In this work we used k — 12 in the Eq. (jSJ). For 
the hep crystal, Q h & cp = 0.484. 



IV. RESULTS FROM SIMULATIONS 



We performed molecular dynamics (MD) simulations in the canonical ensemble using 500 
particles in a cubic box with periodic boundary conditions, interacting with the potential 
Eq. ((H). The parameters employed were a = 5, r /er = 0.7, and c = 1.0. The cutoff radius 
was set to 3.5a^. In order to keep fixed the temperature, the Nose-Hoovep^ thermostat was 
used with the coupling constant qxu = 2. Pressure, temperature, and density are shown in 
dimensionless units, 

Pa 3 

P* = — (6) 



P* = per 3 - (8) 

The translational and orientational order parameters were calculated over 1 000 000 
steps MD simulations, previously equilibrated over 200 000 steps. For low temperatures 
(T* < 0.4), additional simulations were carried out with equilibration over 500 000 steps, 
followed by 2 000 000 steps simulation run. The time step was 0.002 in reduced units. 

For studying the crystalline structure of our model we consider the expected, follow- 
ing conformations for the ramp potential^: simple cubic (sc), body centered cubic (bec), 
face centered cubic (fee), simple hexagonal (sh), hexagonal close packing (hep), and the 
rhombohedral-60° (rh60). Perfect crystals with such conformations were constructed and 
the configurational energy per particle u = U*/N was calculated for each arrangement. In 
the canonical ensemble, the most stable crystal at the ground state is that one with lower u. 
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From the Fig. El we see that the hep conformation is the more stable for densities p* < 0.107 
(see inset). The bee conformation is the more stable one for 0.107 < p* < 0.187 (not shown). 

10 I 1 1 1 1 1 1 r-i 
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FIG. 2: The dimensionless configurational energy per particle for the several crystal structures 
considered: face centered cubic (fee), body centered cubic (bec), simple cubic (sc), simple hexagonal 
(sh), hexagonal closest packing (hep), and rhombohedral-60° (rh60). We see that the hep has the 
lower configuration energy per particle for densities p* < 0.107 (see inset). Hence, the expected 
structure for our model at T* = is the hep for p* < 0.107. For 0.107 < p* < 0.187 the bec phase 
has the lower configurational energy between those studied (not shown). 

Studying the equation of state pressure against density for our model we found a mono- 
tonic behaviour for P(p). Hence, an increase in pressure means increase in density as shown 
in Fig. El 

Results for the translational order parameter for the liquid phase can be seen in the Fig. 
EJ While for a normal liquid t increases under compression, for our system this is the case 
only for high temperatures. For lower temperatures t presents a local maximum at a density 
Pt-max and a local minimum at a density p t -mm > Pt-max for temperatures T* < 1.5. 
Between pt-max and pt-min an unusual behavior for the translational order parameter is 
observed: An increase in density induces a decrease in translational order. This behavior 
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FIG. 3: The reduced pressure as a function of the reduced density. The seven isotherms show that 
the relation between P* and p* is monotonia 

can be understood analyzing the dependence of the radial distribution function (RDF) upon 
density [see Eq. (J2J)]. The (a), (b), and (c) arrows in the Fig. |U correspond to the density 
range spanned by the Fig. Ela) , E[b) , andEfc) respectively. 

The Fig. El shows the RDF for T* = 0.25 and several densities: (a) p* = 0.04, 0.06, 0.07 
and 0.08; (b) p* = 0.1, 0.11, 0.12, 0.14, and 0.16; (c) p* = 0.18, 0.2, 0.22, and 0.24. The 
arrows indicate the directions of increasing p* and the dashed line is the reduced interparticle 
potential shown in the Fig. ^multiplied by a factor of 0.5 just for clarity. From Fig. Efa) 
we see the growth of g(r) at r/a ~ 2.5 upon compression, causing an increasing of t over the 
range 0.04 < p* < 0.08. See the isotherm T* = 0.25 at Fig. 0] In this range of densities the 
particles are repelled by the repulsive shoulder and the most probable separation is about 
2.5 units (corresponding to the edge of the shoulder in the potential). We see that over the 
intermediate density range 0.08 < p* < 0.18, t decreases as the density increases. Looking 
at the Fig. Efb) one can explain why this happens. Both an increase of g(r) at r/a « 1.0, 
approximating the RDF to 1 (what decreases \g(r) — 1|), and a decrease of g(r) in the next 
peak (close to 2.5) upon compression, causes t to decrease. This new peak at about 1 unit 
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corresponds to the position of the hard-core part of the potential^. Finally, t returns to 
increase upon compression for p* > 0.18. The sharp growth of g(r) at rjo ~ 1.0 above 
the unity [see Fig|Sfc)] underlies this behavior, indicating that all the particles are pushed 
together up to their hard-cores. This was the same behavior observed for the RDF of the 
ramp potential^^k Anomalous variations in t are absent for T* > 1.5 because the thermal 
energy washes out the effect of the repulsive shoulder. 




FIG. 4: The translational order parameter t as a function of the density p* . From top to bottom, 
the sixteen isotherms are T* = 0.25, 0.30, 0.35, 0.40, 0.45, 0.50, 0.55, 0.70, 1.0, 1.5, 2.0, 2.5, 3.0, 
3.5, 4.0, and 5.0. The (a), (b), and (c) arrows corresponds to the density range spanned by the 
Fig- Et a ))Eth)) and[5^c) respectively. The bold line indicates the isotherm T* = 1.5. For T* > 1.5 
no anomalous behavior is observed for t. The line connecting the points are just a guide for the 
eyes. 

For a normal liquid, it is expected that the orientational order parameter, Q$, increases 
under compression. For our potential, however, a local maximum is detected for Q 6 at a 
density p Qmax in such a way that p t - max < PQmax < Pt-min ( see figure EJ). This means that 
for densities between pQ max and pt-mm both the structural order parameters t and Qq have 
an anomalous behavior, since t and Qq decrease under increasing of pressure. We call this 
range of densities the structural anomaly domain. 
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x=r/a 

FIG. 5: The radial distribution functions for T* = 0.25 and several densities: (a) p* = 0.04, 0.06, 
0.07, and 0.08; (b) p* = 0.10, 0.11, 0.12, 0.14, and 0.16; (c) p* = 0.18, 0.20, 0.22, and 0.24. The 
arrows indicate the direction of increasing p* . The dashed line is the reduced interparticle potential 
shown in the Fig. Q multiplied by a factor of 0.5 just for clarity. 

The relation between the several anomalies presented for this potential is shown in the Fig. 
|7| The temperature of maximum densities (TMD) and the diffusivity extrema (DE) lines 
were obtained from previous work^. The TMD line indicates the region of thermodynamic 
anomaly region, inside which the density increases when the system is heated at constant 
pressure. The DE lines determinate the region of dynamic anomaly. Inside this region, 
diffusivity increases with increasing density. In this work we determinate additional three 
lines shown at Fig. [7| the curve of t maxima (C), the curve of Qq maxima (B), and 
the curve of t minima (A). We call the region between the curves A and B the structural 
anomalous region, inside which both the order parameters, t and Qq, become anomalous, 
namely, decrease with density. The curve B, composed by the Qq local maxima points, 
terminates at T* « 5.0, not shown in the Fig. [7|for clarity. As the temperature T* tends 
to 5.0, the densities for the Qq maxima loci tends to zero. For T* > 5.0 we have studied 
the temperatures T* = 5.5, 6.0, 6.5, 7.0, and 8.0 (not shown). For all these temperatures 
the same qualitative behavior for Q 6 was observed: The Qq parameter has local maxima at 
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FIG. 6: The orientational order parameter Qq as a function of the density p*. From top to bottom, 
the sixteen isotherms are T* = 0.25, 0.30, 0.35, 0.40, 0.45, 0.50, 0.55, 0.70, 1.0, 1.5, 2.0, 2.5, 3.0, 
3.5, 4.0, and 5.0. For 5.0 < T* < 8.0 (not shown) the Qq local maximum points occur at p* = 
with a global minimum at p* ~ 0.3. We do not study the cases where T* > 8.0 (see the text for 
more details). The line connecting the points are just a guide for the eyes. 

p* = and global minima at p* ~ 0.3. The ratio between these extrema (local maxima and 
global minima) does not extrapolate 3.5% in any case. We do not simulate temperatures 
T* > 8.0. 

For the SPC/E water^, the region of structural anomalies contains (inside) the region 
of dynamic anomalies, and the thermodynamic anomaly region lies inside this last one. For 
the silica, other tetrahedrally bonded molecular liquid, simulations show*^ an inverse order 
between the structural and dynamic anomaly regions: The diffusion anomaly region englobes 
the structural anomaly region that englobes the thermodynamic anomaly region. For our 
model, we see a water-like cascade of anomaly regions similar to that found for the SPC/E 
water (see Fig. |7|). This suggests that the role played by the structure in our potential, like 
in water, is determinant for giving rise to the other anomalies. 

From the Fig. Elits clear that for densities over 0.107 the hep conformation is not the most 
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FIG. 7: The relationship between the several anomalies presented for our model. The curve B 
ends at T* ~ 5.0 and it is not entirely shown for clarity. See the text for more details. Between 
the Qq maxima line (curve B) and t minima line (curve A) both the translational and orientational 
order parameters t and Q$ become anomalous, namely, decrease with density. This region we call 
structural anomaly region. The diffusion extrema (DE) lines enclose the region inside which the 
diffusion decreases with density — the dynamic anomaly region. The temperature of maximum 
density (TMD) line englobes the region that density anomaly appears. Both the DE and TMD 
lines were obtained from previous work^. These cascade of anomalies presents the same hierarchy 
as observed for the SPC/E watei4*^. 

stable one between those studied. The bcc configuration becomes the expected crystalline 
arrangement for densities 0.107 < p* < 0.187. Hence, the calculation of the Q$ parameter 
was repeated using eight first neighbors (k = 8) in the Eq. (j2J). For this new calculation, 
the entire curve B in the Fig. [7|is shifted approximately 13% (not shown) in the direction of 
lower densities, crossing the curve C. In this new scenario, the structural anomaly region now 
lies between the curves C and A, region in which both Q$ and t decrease upon compression. 
Despite of this change, the whole qualitatively result is not modified, once the structural 
anomaly region remains outside the dynamic and thermodynamic regions. 

As discussed in the introduction, the convenient orientational order parameter for tetra- 
hedral liquids^^ is q. It was reported that for SPC/E water the isothermal paths in a t — q 
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FIG. 8: The t — Qq plane, or order map. Each line corresponds to an isotherm and the arrows 
indicate the direction of density growth. From top to bottom, the isotherms showed here are: T* = 
0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.7, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, and 5.0. Unlike the SPC/E 
watei* 8 -, the paths formed by the t and Qq parameters developed a two dimensional region in the 
order map for temperatures and densities inside the structural anomalous region. As observed for 
the SPC/E water£. 

diagram ordermap collapse into a single line in the structural anomaly region 8 -. This prop- 
erty supports the idea that in water the anomalies in translational diffusion and in rotational 
mobility are relate d 3110 ' 11 . 

In order to check if t and Qq are also related in our isotropic model, the order map was 
also constructed. The Fig. |K]shows the behavior of t as a function of Qq. The arrows indicate 
the growth of density for each isotherm. Similar to the results found for the SPC/E water 8 , 
Silica^ 3 -, and for the ramp potentia l 40 ' 41 , was observed an inaccessible region for the order 
map of our model. However, differently from the SPC/E water—, and similarly to the ramp 
potentia l l41 , the parameters t and Qq do not fall into a straight line in the order map for 
densities and temperatures inside the region of structural anomalies (note in the Fig. [HI that 
t and Qq develop a two dimensional region in the order map). 
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V. CONCLUSIONS 



Using molecular dynamic simulations we have studied the structure of fluids interacting 
via a three-dimensional continuous core-softened potential with a continuous force. The 
translational (t) and orientational (Qq) order parameters introduced by Steinhardt el. al. 18 
were analyzed in the framework proposed by Yan el. alM to quantify the structure order for 
an isotropic liquid. 

Our model exhibits a region of density anomaly, inside which the density increases as 
the system is heated at constant pressure, and a region of diffusion anomaly, where the 
diffusivity decreases with increasing density^ 3 *. In the pressure-temperature phase diagram, 
the density anomaly region lies inside the diffusion anomaly one. 

Complementary to the thermodynamic and dynamic anomalies, both t and Qq behave 
anomalously in a large region of the temperature— density plane, as follows. The parameter t 
have both a local maximum, at a density pt-maxi an d a local minimum at a density pt-min > 
Pt-max- For densities in the range pt-max < p < Pt-min the translational order parameter 
decreases under pressure. For normal liquids the opposite behavior is expected. For the 
parameter Qq, a maximum at a density pQ m ax between p t - m ax and pt-min wa s observed. 
Hence, both t and Qq become anomalous for densities in the range PQmax ^ P ^ Pt-min- 

The 

loci of the Qq maxima, t maxima, and t minima were plotted in a temperature— density plane 
and we showed that the region where t and Qq behave anomalously encloses the regions of 
density and diffusion anomalies discussed above. This is the same behavior observed for the 
SPC/E water^. Differently from SPC/E water, the parameters t and Qq do not fall into 
a straight line in the order map for densities and temperatures what suggests that unlike 
water t and Qq are independent in the anomalous region. 

In resume, the studied continuous core-softened pair potential, despite not having long- 
ranged or directional interactions, exhibits thermodynamic, dynamic^ 3 , and structural 
anomalies similar to the ones observed in SPC/E watei"^. Therefore, we can conclude 
that the presence of anisotropy in the interaction potential is not a requirement for the 
presence of thermodynamic, dynamic and structural anomalies. 
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